
options(error=stop)
setwd("~/Desktop/causal res/SNMM_Local/3. Data and Programming")
# rm(list=ls())
library(xtable)
source("liteep.R")

library(abind)
source("2.0.0 MyFunc.R")

# load(paste("./cluster/submitted091516/RDataFiles/results.091516.SEED",1,".RData",sep=""))
# load("./RDataFiles/0820201.RData")

# btstrap.times = 300
# btstrap.each  = 50
# 
{
dim = dim(results);
btstrap.times = dim[1]
# dim[1] = btstrap.times
dimnames = dimnames(results)
dimnames[[1]] = paste("sim=",1:btstrap.times)
dimnames(results)$`Sample size` <- c("n=100","n=200","n=500","n=1000")

Results = array(NA, dim, dimnames)
Results <- results
# Results <- results[1:300,,,,,,]
}
dr.res.10[1,1,,,,"est"]
dr.res.10[,,1,1,"alpha","nll"]
# for(i in 1:10){
#     load(paste("./cluster/submitted091516/RDataFiles/results.091516.SEED",i,".RData",sep=""))
#     btstrap.index = (i-1) * btstrap.each + 1:btstrap.each
#     Results[btstrap.index,,,,,,] = results
# }
# Results <- Results[1:10,,,,,,]
{
means = apply(Results, 2:7, mean)
sds   = apply(Results, 2:7, sd)
mses   = apply(Results, 2:7, sd)

MyAttach(params.true)
rownames(gamma.true) = c("a0","a1")
para.true   = abind(alpha.true, beta.true, along = 3)
dimnames(para.true) = dimnames(results)[4:6]

## comparison of likelihood and computation time
dimnames(results)
means["n=1000",,,,"alpha","est"]
means["n=500",,"1","1","alpha","nll"]
means["n=1000",,"1","1","alpha","time"]
means["n=100",,,,"beta", "est"]
}
{
## table output in the paper
# params.true$contrast.true = GetContrast(data, alpha.true, beta.true) # DEBUG
para.true.comb    = rbind(ftable(para.true, row.vars = c("Param","1-5")),
                          cbind(params.true$contrast.true,rep(NA,3)))
  
comp.mean  = t(means["n=1000",,1:3,1,1,"contrast"]) 


# cols.to.output <- c("PMLE4","DR")
cols.to.output <- c("DR","PMLE4")
# `Sample size`; est.method;`1-5`; `1/2`; Param; Measure
comp.mean  = cbind(comp.mean[,3],rep(NA,3),comp.mean[,2],rep(NA,3))
table.mean = ftable(means["n=1000",cols.to.output,,,,"est"],
                    col.vars = c("est.method","1/2"),row.vars = c("Param","1-5"))
table.bias = rbind(table.mean,comp.mean) - cbind(para.true.comb,para.true.comb)
table.bias = RoundS(table.bias*100,2)
table.sd =  ftable(sds["n=1000",cols.to.output,,,,"est"],
                  col.vars = c("est.method","1/2"),row.vars = c("Param","1-5"))
comp.sd    = t(sds["n=1000",,1:3,1,1,"contrast"])
comp.sd    = cbind(comp.sd[,3],rep(NA,3),comp.sd[,2],rep(NA,3))



# cols.to.output <- c("MLE3", "PMLE4","DR")
# # `Sample size`; est.method;`1-5`; `1/2`; Param; Measure
# comp.mean  = cbind(comp.mean[,1],rep(NA,3),comp.mean[,2],rep(NA,3), comp.mean[,3],rep(NA,3))
# table.mean = ftable(means["n=1000",cols.to.output,,,,"est"],
#                     col.vars = c("est.method","1/2"),row.vars = c("Param","1-5"))
# table.bias = rbind(table.mean,comp.mean) - cbind(para.true.comb,para.true.comb, para.true.comb)
# table.bias = RoundS(table.bias*100,2)
# table.sd =  ftable(sds["n=1000",cols.to.output,,,,"est"],
#                    col.vars = c("est.method","1/2"),row.vars = c("Param","1-5"))
# comp.sd    = t(sds["n=1000",,1:3,1,1,"contrast"])
# comp.sd    = cbind(comp.sd[,1],rep(NA,3),comp.sd[,2],rep(NA,3),comp.sd[,3],rep(NA,3))
# 
# 
table.se = RoundS(rbind(table.sd,comp.sd)*1000/sqrt(btstrap.times),2)
table.bias.se = NULL
for(j in 1:ncol(table.bias)){
    table.bias.se = cbind(table.bias.se,
                          paste(table.bias[,j],"(",table.se[,j],")",sep=""))
}
table.bias.se[11:13,c(2,4)] = "$-$"
table.bias.se = table.bias.se[c(1:5,7,6,8,10,9,11:13),]

# row.names = c(paste0("$\\", c("theta_1$","theta_2$","theta_3$","theta_4$","theta_5$",
#               "phi_0$","phi_1$","gop$","eta_0$","eta_1$")),
#               c("$E[Y(0,0)]$",
#                  "$E[Y(1,1)]$", "$E[Y(1,1)]/E[Y(0,0)]$") )

row.names = c(paste0("$\\", c("theta(\\ast)$","theta(1,1)$","theta(1,0)$","theta(0,1)$","theta(0,0)$",
                              "phi_0$","phi_1$","text{gop}$","eta_0$","eta_1$")),
              c("$E[Y(0,0)]$",
                "$E[Y(1,1)]$", "$E[Y(1,1)]/E[Y(0,0)]$") )

 # c("theta_1$","$theta_2$","$theta_3$","$theta_4$","$theta_5$",
 #              "$phi_0$","$phi_1$","$gop$","$eta_0$","$eta_1$","$E[Y(0,0)]$",
 #              "$E[Y(1,1)]$", "$E[Y(1,1)]/E[Y(0,0)]$") 
rownames(table.bias.se) = row.names
# table.bias.se
# xtable(table.bias.se)
print(xtable(table.bias.se), sanitize.rownames.function = identity)
# save(Results, means, para.true, table.bias.se, file="./RDataFiles/results091516.RData")
}

#####################################################################





t.misgop <- table.bias.se
t.logit <- table.bias.se


print(xtable(cbind(t.misgop[1:5,1:2], t.logit[1:5,1:2])), sanitize.rownames.function = identity)
print(xtable(cbind(t.misgop[1:5,3:4], t.logit[1:5,3:4])), sanitize.rownames.function = identity)


# Print the table of computation time
{
  cluster.path = "~/Desktop/causal res/SNMM_Local/3. Data and Programming/cluster/080820"
  setwd(cluster.path)
  load(paste("./RDataFiles/16Jul2021-n-",1000,"-theta-not-0-l0-cont",1,".RData",sep=""))
  btstrap.times = 300
  
  dimnames(results)$`Sample size` <- c("n=100","n=200","n=500","n=1000")
  # Results <- results
  Results <- results[1:300,,,,,,]
  ## table output in the paper
  # params.true$contrast.true = GetContrast(data, alpha.true, beta.true) # DEBUG
  means = apply(Results, 2:7, mean)
  sds   = apply(Results, 2:7, sd)
  cols.to.output <- c("MLE3", "PMLE4","DR")
  time.mean.cont  = t(means["n=1000",,1,1,1,"time"]) 
  
  load(paste("./RDataFiles/23Jul2021-n-1000-theta-not-0-l0-disc",1,".RData",sep=""))
  Results <- results
  means = apply(Results, 2:7, mean)
  sds   = apply(Results, 2:7, sd)
  time.mean.discrete  = t(means["n=1000",,1,1,1,"time"]) 
  # row.names = c(paste0("$\\", c("theta_1$","theta_2$","theta_3$","theta_4$","theta_5$",
  #               "phi_0$","phi_1$","gop$","eta_0$","eta_1$")),
  #               c("$E[Y(0,0)]$",
  #                  "$E[Y(1,1)]$", "$E[Y(1,1)]/E[Y(0,0)]$") )
  
  row.names = c("$L_0 \\text{discrete}$","$L_0 \\text{continuous}$")
  table.time <- round(rbind(time.mean.discrete, time.mean.cont),2)
  rownames(table.time) <- row.names
  # c("theta_1$","$theta_2$","$theta_3$","$theta_4$","$theta_5$",
  #              "$phi_0$","$phi_1$","$gop$","$eta_0$","$eta_1$","$E[Y(0,0)]$",
  #              "$E[Y(1,1)]$", "$E[Y(1,1)]/E[Y(0,0)]$") 
  # rownames(table.bias.se) = row.names
  # table.bias.se
  # xtable(table.bias.se)
  print(xtable(table.time), sanitize.rownames.function = identity)
  # save(Results, means, para.true, table.bias.se, file="./RDataFiles/results091516.RData")
}



################
###########
#######
{
  nrow(marginal.data)
  idx.00 <- which(marginal.data$a0a1=="00")
  idx.01 <- which(marginal.data$a0a1=="01")
  idx.10 <- which(marginal.data$a0a1=="10")
  idx.11 <- which(marginal.data$a0a1=="11")
  g.val.dat <- cbind(marginal.data$value[idx.00], marginal.data$value[idx.01],
                     marginal.data$value[idx.10],marginal.data$value[idx.11])
  colnames(g.val.dat) <- c("00","01","10","11") # two numbers are (a0,a1)
  to.output <- cbind(log(g.val.dat[,"10"]/g.val.dat[,"00"]),
                     log(g.val.dat[,"11"]/g.val.dat[,"10"]),
                     log(g.val.dat[,"01"]/g.val.dat[,"00"]),
                     log(g.val.dat[,"11"]/g.val.dat[,"00"]))
  means = round(apply(to.output, 2, mean),4) * 100
  se = round(apply(to.output, 2, sd)/ sqrt(500),4)* 100
  
  # to.output <- cbind(g.val.dat[,"10"]-g.val.dat[,"00"],
  #                    (g.val.dat[,"11"]-g.val.dat[,"10"]),
  #                    (g.val.dat[,"01"]-g.val.dat[,"00"]),
  #                    (g.val.dat[,"11"]-g.val.dat[,"00"]))
  # means = round(apply(to.output, 2, mean),3)
  # se = round(apply(to.output, 2, sd)/ sqrt(500),3)
  
  table.bias.se = NULL
  for(j in 1:length(means)){
    table.bias.se = cbind(table.bias.se,
                          paste(means[j],"(",se[j],")",sep=""))
  }
  table.bias.se = t(table.bias.se) 
  row.names= c(paste0("$\\", c("log E[Y(1,0)]/E[Y(0,0)]$",
                                "log E[Y(1,1)]/E[Y(1,0)]$",
                                "log E[Y(0,1)]/E[Y(0,0)]$",
                                "log E[Y(1,1)]/E[Y(0,0)]$")))
  
  # c("theta_1$","$theta_2$","$theta_3$","$theta_4$","$theta_5$",
  #              "$phi_0$","$phi_1$","$gop$","$eta_0$","$eta_1$","$E[Y(0,0)]$",
  #              "$E[Y(1,1)]$", "$E[Y(1,1)]/E[Y(0,0)]$") 
  rownames(table.bias.se) = row.names
  # table.bias.se
  # xtable(table.bias.se)
  print(xtable(table.bias.se), sanitize.rownames.function = identity)
}




